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Abstract 

We propose an algorithm for the orthogonal fast discrete spherical Bessel 
transform on an nniform grid. Onr approach is based npon the spherical 
Bessel transform factorization into the two snbseqnent orthogonal transforms, 
namely the fast Fonrier transform and the orthogonal transform fonnded on 
the derivatives of the discrete Legendre orthogonal polynomials. The method 
ntility is illnstrated by its implementation for the nnmerical solntion of the 
three-dimensional time-dependent Schrodinger eqnation. 
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1. Introduction 


The discrete spherical Bessel transform (DSBT) arises in a nnmber of 
plications, snch as, e.g., the analysis of the cosmic microwave backgronnd [l|, 
the nnmerical solntion of the differential eqnations 2|, y, 1^, and the nnmeri- 
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cal evaluation of multi-center int^rals { 5 , 6 ]. Many different SBT algorithms 


have been proposed so far 


t 


10l |. But none of them possess all of the 


advantages of their trigonometric progenitor, namely the fast Fourier trans¬ 
form (FFT). These advantages are the performance fastness, the uniform 
coordinate grid, and the orthogonality. 

An example of the problem requiring the simultaneous presence of all the 
advantages is the solving of the Schrodinger-type equation (SE) by means 
of the pseudospectral approach [9j. The grid uniformity provides the same 
accuracy of the wave function description in the whole domain of dehnition. 
The grid identity for all orders of a spherical Bessel functions (SBF) allows 


to switch to the discrete variable representation (DVR) 


The DSBT or¬ 


thogonality is needed to provide the hermiticity of the radial part of the 
Laplacian operator in DVR. The lack of the Laplacian operator hermiticity 
impedes the convergence of iterative methods (such as conjugate gradient 
method) for the solution of matrix equations (which are obtained by DVR 
from the stationary SE). In the case of time-dependent SE, the hermiticity 
of the Laplacian operator is crucial for the conservation of the wave function 
norm during the time evolution. 

A pioneering approach based upon the convolution integral 
requires a number of operations of the order of N log 2 N for its performing, 
just like the FFT does, that means that it is quite fast. However it em¬ 
ploys a strongly nonuniform grid (a node location exponentially depending 
on its number). Hence the attempts of its utilization for the SE solving 
2|, l3| ended in problems with the strong near-center localization of a wave 
function. A method rest on the spherical Bessel functions expansion over 
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the trigonometric functions js] also appears to be quite fast (requiring as 
few as {i + l)iVlog 2 operations) and employs a uniform grid. But it is 
not orthogonal and has stability difficulties because of the singular factors 


in the spherical Bessel functions expansion over the trigonometric 
Next, a Gauss-Bessel quadrature based technique suggested in jo, 


unctions. 


13| is or¬ 


thogonal and converges exponentially, but it is not fast (as the number of 
operations required scales as iV^) and needs an ^dependent grid. Neverthe¬ 
less its fast convergence and the near-uniform grid motivated to apply it for 


a time-dependent Gross-Pitaevsky equation 


Finally, an approach rest on 


the SBF integral representation via Legendre polynomials, proposed in 


io|, 


appears to be fast, makes use of the uniform grid, but it is not orthogonal. 

In the present work we are proposing the algorithm for the DSBT that 
is orthogonal, fast, and it implies the uniform grid. Our approach is based 
upon the SBT factorization into the two subsequent transforms, namely the 
FFT and the discrete orthogonal Legendre polynomials derivatives based 
transform. 

The paper is organized as follows. In Section [2l we develop the orthogo¬ 
nal fast DSBT on a uniform grid. Next, in Section [3l the proposed method 
is tested via the evaluation of the Gaussian atomic functions transform and 
also the DSBT basis functions comparison to the exact SBFs. In Section 0] 
the DSBT- and DVR-based approach (DSBT-DVR) for the time-dependent 
Schrodinger equation (TDSE) solving is suggested and examined. The ap¬ 
proach efficiency is illustrated by treating of the problem of the Hydrogen 
molecular ion ionization by laser pulse. Finally, in Section |5] we briefly dis¬ 
cuss the obtained results as well as the prospects of DSBT and DSBT-DVR 
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application. 


2. Development of the method 
2.1. Basic formulation 

A typical problem involving the spherical Bessel transform (SBT) is the 
plane wave expansion of a three-dimensional function \&(r, 6, 0). The expan¬ 
sion over the spherical harmonics yields a radius-dependent function 

= j 0 )T(r, 6, ct>)dn. (1) 

If the function T(r, 0, 0) has no singularities, then —)■ 0) ~ r^. 

Let us introduce the SBT as 

= \l^ xe{kr)ife{r)dr (2) 

Here we perform the function substitution 0£(r) = (a magnetic 

quantum number is not used further, therefore from this point on we omit it 
from the denotation for the sake of simplihcation), then execute the expansion 
over the functions 


Xe{x) = xji(x), 


( 3 ) 


where je{x) is a spherical Bessel function (SBF) of the hrst kind. The 
functions Xi{kr) satisfy the normalization condition Xi{^'^)xe{^''^)dr = 
{7i/2)6kk'- The pre-integral factor in ([2]) is introduced in order to make the 
transform ([2]) unitary. 

The beginning of our derivation coincides with the one in the work 10 |. 


But unlike its authors we are going to aim at the factorization of the SBT into 
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the two separate transformations, namely the FFT and also the additional 
orthogonal transform which we denote Fourie-to-Bessel transform (FtB). The 
SBF may be presented as 


je{z) = ^ I Piiv) exp{iz7])dT] 


(4) 




where Pe{f]) is the Legendre polynomial of f'-th order. Upon substituting the 
latter expression into Eq.([2]), we obtain 


ce{k) = J' Piiv) krP^^^Mr)drd7i 

Here the integral over r is different from the Fourier transform of the func¬ 
tion by the presence of the integrand factor kr. This factor might be 

represented as a result of taking a derivative of over r]. Thus we get the 
expression 


/T 1 /■! 8 r°° 

Making use of the Legendre polynomials parity condition = {—lyPiiri), 

one may further reduce the integral over rj from —1 to 1 to the one in the 
limits from 0 to 1 as 


= Vf^ r^''" ■ (-P~''‘'''']Mr)drdv (5) 

Next, let us dehne a new function 

em = yf ^ jyp' - (-l)'e-“1*(r)ir (6) 

The term — (—l)^e“®^®']/(2i^+^) is equal to sin(fcr) for the even 

i and to (—l)F/2l cos(fcr) for the odd ones. Hence the expression (E]) appears 
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to be correspondingly the sine/cosine Fourier transform of the function 'ipe{r), 
depending on i being even or odd. 

In terms of the new denotation the formula ([5]) can be rewritten as 

ce(k) = 

Upon making the substitution rj = q/k this expression takes the following 
form 

ce{k) = [ Pi>{q/k) ^^^j^^ dq 

Jo dq 

Let us perform the integration by parts, then move the derivative over q to 
the Legendre polynomial. As a result we obtain the following formula for the 
FtB 


ce{k) = ci{k) - h{q)dq (7) 

It is easily seen that c^ik) = Co(fc), just as expected, since Xoikr) = sm{kr) 
and the Bessel expansion coincides with the Fourier expansion at £ = 0. 

One may rewrite d?]) in the operator form as q(A;) = Tce{k), where the 
integral transform operator T has the kernel 

T(k,q) = S{q-k)-e{k-q)Sk^. ( 8 ) 


Here 


6{x) 


0, a: < 0; 

< 1/2, X = 0; 

1, X > 0. 

V 


(9) 
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is the Heaviside step function. The FtB operator T must be unitary (that is 
T~^ = T^), hence the inverse transform is C(,{k) = T'^Ci{k). The inverse FtB 
might be explicitly dehned as 


ci{k) = ce{k) - 


P'ik/q) 


ce{q)dq. 


( 10 ) 


The snbstitution of ffTOjl into ([7]) demonstrates that ffTOj) is indeed inverse 
in respect to ([7]), dne to the condition 


ki 


-dq = 


_ fe) + h(hZM«(fe _ k,). 

ki fc2 


( 11 ) 


This condition holds trne since for any polynomial p{x) of the order s < i 
true is the expression 


P^{x)p{x)dx = p{l) — {—1 Yp{—1) 


( 12 ) 


'-1 


In turn, this relation is the consequence of the well-known Legendre polyno¬ 
mials property, 

ri 


PYx)x^dx = 0 ; p < 


(13) 


'-1 


following from their orthogonality. 


2.2. Discretization of the transform 

Let us introduce the coordinate grid with the step Ar in the following 
way 


p = (i — l/2)Ar; i = 1,..., N. 


(14) 
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The elements of vector ij) of the function values are sampled on the grid as 


i’i ^ ipe(ri)VAr. 


(15) 


The elements of the Fourier transform matrix are dehned as follows 

^ ^iknr _ 

1 

v/(2 - 5no)N 

Here the momentum grid is 

kn = nAk] n = Pi,. .. ,Ne, (17) 

where the momentum step is Ak = Tr/rmax, the integration interval size is 
?"max = NAr, and the summation limits are = [1 + (— 1)^]/2 and Ni = 
N + Pi — I, that is n = 0, ..., iV — 1 for the odd i and n = 1, ..., iV for the 
even ones. 

The Fourier transform may be written in the matrix form as 

f = FV>. (18) 

It yields as a result the value of the vector of the Fourier expansion coefficients 
f related to the function Ci{k) as follows: 

fn Ci(^kn)y^Wn ( 19 ) 


X 


sin(/c„rj), even i; 
cos^knTi), odd i. 


(16) 


where the weights 



Since the Fourier transform matrix F is orthogonal, then under the transform 
the norm is conserved, that is f^f = The transform fflSj) performing 

through the FFT algorithm requires the number of operations of the order 
oiNlog^N. 

The transform ([7]) conserves the norm according to 


poo poo 

/ \ce{k)\‘^dk = / \ce{k)\‘^dk 

'o Jo 


( 21 ) 


The approximation of the integrals in this relation by the trapezoidal rule 
yields 


b^b = ftf. 


( 22 ) 


where we introduce the vector b composed of the coefficients of the Bessel 
expansion 


bn ( kn ) \/ Wn- 


(2.3) 


Next, let us write the direct and inverse discrete FtB (DFtB) in the 
following form 


b = Tf 
f = T ^b. 


(24) 

(25) 


In order for 022 p to hold true, the matrix T has to be orthogonal, that is 

T ^ = T^. (26) 

If we attempt to apply the trapezoidal rule directly to ([7]), then we would 
obtain 

rp _ X Pi{km/kn) 

J-nm — ^nm j \/^ 

Kin 
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However this technique of the matrix construction does not provide its or¬ 
thogonality. The reason is that the equation flT^ does not hold upon the 
approximate integration. The employing of the high-order Newton-Cotes 
rules instead of the trapezoidal rule does not make the situation better. In 
order for flTT]) to be true, it is necessary for flT^ to hold for all the subgrids 
with the arbitrary nodes number. The high-order Newton-Cotes rules do not 
provide high accuracy for an arbitrary subgrid. Therefore the only way to 
preserve the transform orthogonality appears to be the modihcation of the 
integral ([7]) kernel under the proceeding to the numerical integration. 


2.3. Discrete Legendre orthogonal polynomials 

In the context of the summation on grids, the properties analogous to 
those of the Legendre polynomials are possessed by the so-called discrete 
Legendre orthogonal polynomials (DLOP) Q. DLOP Peii,N) satisfy the 


orthogonality property given by 


N 




(27) 


i=0 


and also the normalizing condition Pi{0,N) = 1. Here 

Af(F N) = (^ + ^ +1)— 798! 

where is j-th falling factorial of i, = i{i — 1)... {i — j + 1). DLOP might 
be presented as 

e 

P,(i,N) = (29) 


j=0 


where l{i,j) are coefficients of the expansion of the shifted Legendre polyno¬ 
mial Pi{l — 2x) = ■ This means that Pi{i,N) can be obtained 
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from P^{1 — 2x) through the substitution of ii/Nl for xK As the grid size 
increases, DLOP tend to the usual Legendre polynomials according to 

Pe{t,N) = P,{l-2t/N) + 0{N-^). (30) 

Due to the orthogonality condition fl27|) DLOP possess a property similar 
to the property flT^ . as follows 

N 

'^Pt{i,N)P = 0; s<^. (31) 

i=0 

Making use of this fact one can easily prove (as shown in the Appendix) 
that for any discrete polynomial p{i) of the order fi < i true is the following 
relation 

N 

Y,P'S.N)p{i)w,{N) = (-l)V(iV)-p(0) (32) 

i=0 

where the weight function coincides with the weights of the trapezoidal inte¬ 
gration rule for the grid with a unit step 

Wi{N) = 1-^-. (33) 

Here we dehne a new discrete polynomial 

(34) 

which is proportional to the backward difference 

V[P,](*, iV - 1) = P,(z, iV - 1) - P,(z - 1, iV - 1) (35) 

and hence has the order ^ — This polynomial tends to the derivative of the 
usual Legendre polynomial as 

P;(*,7V) = lp,(i-2z/iV) + 0(iV-3), (36) 
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that is faster than the DLOP tends to its non-discrete analogue. Therefore 
we shall further refer to P^{i, N) as the derivative of the discrete orthogonal 
Legendre polynomial (DDLOP). It should be mentioned that this term has 
different meanings throughout the literature jl^ . 

It follows from the relation fl5^ that the integral kernel in Eq.(I7]) can be 
approximated on the grid ffT7|) by means of DDLOP according to 

Peik ^ ^ 0[{Ak)-% (37) 

2.4- Transform matrix 

Let us suppose the transform matrix T has the elements as follows: 


Tnm Oln [<5nm Lnm\ j 


(38) 


where the lower triangular matrix L is an approximation of the kernel of the 
integral (j7]) via ([37]), dehned by 

Tnm = - m)P^{n - m, 2n)\jl- (39) 


Here the Heaviside function 6{n — m) is specihed in (I9]), and the additional 
factor \J^ ~ provides the weight function fl20l) in the product Tf . Dehned 
in such a way Tnm exist only for n > hqi where 


noi — 


'i + 1 
2 


(40) 


since there are no DLOP of the order i at smaller n. 

By making use of Eq. fl32l) . one may show (see the Appendix for details) 
that the rows of the matrix I — L are mutually orthogonal, that is 

Nt 

^ ^ [^nZ [5,7^; Lmll ^nm-, ^ ^ '^01 

l=Pt 
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The rows of the matrix T are equal to those of the matrix I — L multiplied 
by the normalizing constants 



- 1/2 


(42) 


Here the normalizing constants 1 at n —)• cxo. 

Since the approximation flTTIl implies the error scaling as 0[{Ak)^], the 
result of the DFtB defined by matrix fl38|l differs from the result of the exact 
FtB by the value 0[(AA;)^]. 

2.5. Completion of the basis 

The Eqs. fl38ll39p define only N — hqi + 1 rows in the transform matrix. In 
order to make the basis complete we have to supplement it by extra hqi — 1 
vectors that are orthogonal to all other ones. 

The DDLOP property Eq. fl3^ (which Eq. fITT]) follows from) leads to the 
fact that the basis vectors specified via Eqs. fl38ll39|) are to be orthogonal to 
any polynomial of the order s < l — \. That is to say, one might construct the 
extra basis vectors from the polynomials of the order s < i — 1. To provide 
these extra vectors to be orthogonal to not only the basic basis vectors, but 
to each other as well, we shall choose them as the DLOPs (not the DDLOP!) 
of the corresponding orders. 

So we shall define the extra basis vectors as follows 


nm 



(43) 


where the polynomial order is 


l{n) = 2n-pi, 


(44) 
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and the normalizing constant is where M is given in 

Eq.(l28]). 

Now let ns elncidate the physical meaning of the extra basis vectors (jl3]). 
For this pnrpose we shall derive their appearance in coordinate representa¬ 
tion. First, let ns begin with the expression for the coefficient 

I T 

bn = (45) 

m=pi 

and proceed making nse of the relation between the FLOP and the nsnal 
Legendre polynomial 

Pi{N — m,2N) ^ Piim/N). 

By changing the snmmation to integration we get 

P^max 

bn ~ «n / Pi(n)iq/kN)ceiq)dq, 

Jo 

where = NAk. Finally, snbstitnte here (E]) and nse Fq.(jl]) to obtain 

bn ~ I JKn){kNr)Mr)dr (46) 

Thns, bn Sit n < nor are the coefficients of the expansion in the fnnctions 
3 i(n){kiqr) (while at n > nor are the coefficients of the '0r(’") expansion in 
Xeiknr))- At r —)■ 0 asymptotics are Xi{kr) n-u and ji(n){kNr) r'-JJ. 
Since /(n) < £ — 1, the additional basis vectors represent the high-momentnm 
wave-fnnction components converging to zero slower than That is to 

say, the additional vectors emerge to be a linear combination of the non- 
regnlar SBFs (SBFs of the second kind). 

If '0r(r) appears to be a resnlt of the spherical harmonics expansion of 
a 3D fnnction having the continnons derivatives np to order ^ + 1, then 
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i^e{r 0) hence it should be = 0 at n < hqi. In the case when 

the function possesses singularities such components become non-zero and 
should be considered as having the energy larger than that of the component 
specihed by the coefficient b^- 

Thus, our DSBT does not yield coefficients for the small momenta kn, n < 
noi- The reason for this is that at such kn there exist no SBF of the first 
kind satisfying the boundary conditions (see also the Section |3]) . 


2.6. Fast multiplication by transform matrix 

The fast transform might be accomplished by means of the technique 
proposed in the work 1^, except that we shall use DDLOP instead of the 


Legendre polynomials. 

First, let us expand DDLOP over the powers of the variable m in the 
following way: 

i-i 


P^{n — m, 2n) = — ^ 


u=0 


Upon substituting this expansion into fl2T|) one obtains 


e-i 


bn (^nfn Q^n ^ ^ 


u=0 


Here we introduce the notation 

5: 


n—1 


^UTl. 


^/o + rn^fm + - 

^ m=l 


rFfn. 


The sum flT^ can be evaluated through the recurrence relation 


(47) 


(48) 


(49) 


^un ^u,n—l T 2 fn T ('^ 1) fn—l\ 


(50) 
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Thereby the sums calculation for all the n G [uq^, N] requires as less as N 
operations. In total one needs to evaluate |'^/2] x (iV + 1 — riQi) sums (|'^/2] 
appearing because of the fact that the half of the coefficients in flTTl) are zero 
due to the DDLOP parity), then to perform the summation over u in flTHj) . 
resulting in the altogether operations number scaling as iN. 

At n < nog the coefficients bn are to be computed according to fHSD . That 
means that for every n < nog one has to calculate the vectors scalar product 
that requires extra 0{£N) operations, that is the operations number scaling 
is the same as for the bn evaluating for n > nog via fl50ll48p . In sum, to 
accomplish the transform fl241) one needs to perform 0{iN) operations. 

Now let us consider the inverse transform. The substitution of flTT)) into 
fl2S]l yields for m > nog the following: 

noe-l 1-1 

fm ^ ^ Tnmbn T Oimbm ^ ^ ^ ^urm (^1) 

n=p£ u=0 

where 

Ni ^ 

^um ^ ^ Otn^guij^^bn T '^Cini^guijA^bm- (^2) 

n=m+l 

Next, for m < nog we obtain 


fm = 


noe-l 

n=pt 


1-1 


v=0 




-Tfl + 


(63) 


where 


The sum (|52|1 may be evaluated according to the recurrence relation 

^um ^u,7n-\-l T 2 fgv{m)bm + am+iigu{m + l)bm+i] . (55) 
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Thus, the inverse transform 0251) performing requires the number of opera¬ 
tions 0(£iV), just as the direct transform does. 

The fast Fourier transform fflSl) operations number scales as N log 2 N. 
That is the DSBT in total requires 0(iVlog 2 A^) -|- 0{£N) operations. At 
large N the FFT strongly over-demands the DFtB, hence the overall DSBT 
algorithm processing time appears to be defined by the FFT processing time. 

3. Numerical test of the method 

The method convergence was examined on three grids 014p with the same 
space step Ar = 0.4 and the various values of the space region rmax = 
NAr = 51.2, 102.4, and 204.8. That is, the grids possessed the same maximal 
momentum = AkN = vr/Ar and various momentum steps Ak = vr/r ma x. 

Let us begin with the check of the convergence of our transform for 
smooth functions, which are commonly used in atomic physics namely Gaus¬ 
sian atomic orbital functions 

'ip£{r) = exp(—r^/2). (56) 

The FigHlshows the absolute value 6 = |c„£ —Q(fc„)| of the difference between 
the exact SBT result C£{k) (obtained by the numerical evaluation of the 
integral in Eq.(|2])) and the DSBT result Cn£ = it is seen that 

upon the r^ax increasing the error decreases as 1/r^ax (or, equivalently, as 
Afc^). It coincides with the convergence rate expected from theory. 

In order to reveal the main error source, let us perform the compar¬ 
ison of the exact SBFs with the DSBT basis functions in the coordinate 
representation, which correspond to the transform basis vectors. For this 
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Figure 1: Deviation of the Gaussian atomic orbitals expansion from exact for rniax=51.2 
(solid lines), rniax=102.4 (dashed lines) and rniax=204.8 (dotted lines) 

purpose we apply the inverse transform to the vector composed of coeffi¬ 
cients bm = {rmax/‘^Y^‘^Wn^^'^6mn{k), where n{k) = k/Ak and k is the hxed 
momentum. The result function may be written as follows: 

Xm{n) = (57) 

The factor (r ma x/2)^/^ is introduced in order to provide the convergence 
XntiTi) Xe{knri) at r^ax —>■ oo. The results are presented at the FigJU 
For every i we chose the momentum k in such a way that for the maximal 
Ak grid a basis vector number n{k) = The north function has the nodes 
number minimal among the functions satisfying the boundary conditions at 
r = Tmax and the asymptotics ~ at r = 0. Upon the step Ak decreasing 
at the hxed k the number n{k) grows and the DSBT basis function in its 
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1.5 n 


Nl; k=0.061359 


1.5 n 


1=2; k=0.122718 





Figure 2: Convergence of the DSBT basis functions to SBFs with decreasing of Ak. 


domain of definition becomes closer to the SBF. 

The Fig J3] demonstrates the DSBT basis functions at n < riQi. One can see 
these functions behavior is just as expected from (06]), that is they are high- 
frequency, decrease with distance and have asymptotic behavior different 
from ~ at r = 0. 

In order to demonstrate the main error source, we plot the dependence 
of 5 = \Xn(k}i{'>^) — \ oil radius (FigJl|). At that we take the momentum 

value such that n{k) = 9 even at the smallest grid. It is apparent that the 
difference grows with distance. The reason for the difference to be maximal at 
r = risf may be understood if one recalls that SBFs asymptotically equivalent 
to 

Xe{kr) sm[kr — £7r/2 — l{l + l)/2kr]] r —)■ oo. (58) 
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1=3 



Figure 3: Extra basis functions for ^ = 3. 

The DSBT basis function satisfies the same boundary conditions as the func¬ 
tions constituting it (sines or cosines) do 

Xn(fc)f(r) ~ sin(/cr - £7r/2); r-)■ tat. (59) 

That means that there does exist a phase shift /(/ -|- l)/2/crjv near the outer 
boundary between the approximate and exact functions. 

The fact of the amplitude of the difference between the SBFs and DSBT 
basis functions growing roughly linearly with r increasing (which is seen 
from FigH]) indicates that the phase shift between these functions grows 
linearly with r as well. Such phase shift behavior may be interpreted as 
a consequence of the approximate and exact functions wavenumbers being 
distinct. Therefore the DSBT basis function has to be closer to the SBF 
at such momentum k^i. that provides the exact SBF coincidence with the 
DSBT basis function on the grid boundary. This condition might be written 
mathematically as 

Xi{knirmax) = 0, cven i; 

(60) 

Xn^knl'^max') 0, odd £. 
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Figure 4: Difference between the basis functions and SBFs 6 = \Xn(k)e{''') ~ for the 

fixed k = 0.490873843. 


From the phase shift between the DSBT basis functions and exact SBFs near 
the outer boundary one can obtain the approximate expression 

kni ^kn- ( 61 ) 

On FigJSl we plot the 6 = \Xn{k)e{r) — Xe{kner)\, i.e. the difference between 
the DSBT basis functions and SBFs for the corrected momentum kn£- It is 
seen that the difference is an order of magnitude less than in the case of 
the non-corrected momentum (FigSj) at the same r^ax- Yet the rate of 
the basis function convergence to the corrected momentum SBF is still 
(whereas r <C Tmax)- Thus our transformation result arrears to be more 
exact approximation to the expansion in SBFs on the grid kni than it is on 
the uniform grid kn- However the grid kni is non-uniform and ^-dependent, 
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Figure 5: Difference between the basis functions and SBFs S = \Xn{k)e{'i~) ~ Xiikne.r)\ for 
the fixed k = 0.490873843. 

and this complicates the application with no convergence rate advantage. 


4. An example: Solution of the 3D time-dependent Schrddinger 
equation 

4 . 1 . Numerical method 

As an example of DSBT application we developed the method for a nn- 
merical solntion of 3D time-dependent Schrddinger eqnation (TDSE). 

Let ns make nse of discrete variable representation (DVR). We shall begin 
with the fnnction discretization on the 3D grid in the spherical coordinate 
system 

i>ijk = arccosr]j, (j)k)riy/ArAr]jA(j) (62) 
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Here 0^ = A0(fc — 1), A; = 1,..., is a polar angle grid; A0 = 2tt/N^ is a 
grid step; rjj and Ar/j, j = 1,..., are Ganss-Legendre qnadratnre nodes 
and weights correspondingly. Upon DVR implementation, the transforma¬ 
tion given by Eqs. ffT]|^ is written as 

c = BY-j/j. (63) 

Here c is a vector of coefficients of the expansion in spherical waves 

^nlm Qm(A^n) (6^) 

Since fnrther we will need identical momentnm grids for all the angnlar mo¬ 
menta whereas the grid kn dehned by Eq. (lT7|) differs for odd and even i, 
we shall now introduce a new grid for the momentum radial component 

kn = nAk] n = 1,..., A, (65) 

and also shall assume cim{kN)y/wN = cozm for even 1. Validity of this proce¬ 
dure may be justihed as follows. Due to the subsection 12.51 the coefficient 
having n = 0 for all I corresponds to a projection onto a non-regular high- 
frequency function. In turn, a non-regular high-frequency function may be 
approximated by the sum of regular functions with large momenta. So, as the 
highest spectrum part coefficients are evaluated rather inaccurately, one may 
assume without loss of accuracy the coefficient with n = 0 to be the value of 
the projection on a basis function with large momentum kN- Meanwhile, for 
a smooth \b(r, 6 , 0) this coefficient vanishes anyway. 

The rest of designations used in Eq. (l6^ are as follows: B is the SET 
matrix with elements 

BnlmiVm' l\nAllAmm' (66) 
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and Y is the matrix of the transform to the expansion in spherical harmonics 




(67) 


which might be represented as the product of the matrix with elements 


^jmj'k ^jj' 


^i7TL(p}~ 




\/~K^ 


( 68 ) 


and the matrix 


Plmjm' Pi (^j ) \/ 


(69) 


Here P; [r]) are associated 


Legendre polynomials ortlionormalized on the 


15]. The transform Eq. fl6^ requires O{NrN 0 N'^) + 


Gauss-Legendre quadrature 
0{NrN^Ng) + O{N 0 Nfi,Nrlog 2 Nr) operations. It should be noted that the 
term 0{NrNff,Ng) is caused not only by the operation of polar angle P mul¬ 
tiplication by the transform matrix, but also by the multiplication by ma¬ 
trices T; (matrix multiplication at the hxed I requires 0 {NrN^l) operations, 
whereas the number of different Vs is equal to N 0 ). Therefore if one had tried 
using (instead of DVR) any methods that do not employ P transformation, 
it would not make sense, because it would not imply getting rid of the oper¬ 
ations number quadratic growth with N 0 increasing. However, the transform 
algorithm is easily parallelizable, hence may be run in quite modest amount 
of computer time even for the large value of N 0 . 

Besides, we introduce the matrix 


^nlmn'jk ^nn'^ 


(70) 


of the transform to the expansion in terms of modihed spherical harmonics 
|i^ related to the common spherical harmonics as Yim{ 0 ,(j)) = VYimid,(l)) {i 
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is the imaginary unit here). It can be used to perform the switch to the 
momentum DVR 

= Y^BYi/j, (71) 

where the vector ip components relate to the plane wave expansion coeffi¬ 
cients as follows 

'-Pnjk = arccos r]j, 0fc)y/AfcApjA0. (72) 

The very possibility of the transition to the plane wave expansion is the 
main reason for us having introduced the unihed grid for the momentum 
radial components, Eq. fl65|) . 

One may write the Hamiltonian for a particle in external held as 0 

H = Y - + (73) 

Here p = —-iV is the momentum operator, U{r,6,(j),t) is the potential of 
electron-nuclei interaction, and A{t) is the vector potential of the external 
electric held. The vector potential dehnition 

A{t) = — [ qS{t')dt, (74) 

Jo 

slightly dihering from the commonly used one, will be used further for the 
sake of expressions brevity. Here g is a particle charge and S {t) is the external 
electric held strength. 

Employing DVR and the transform Eq. (I63|) allows to represent the Hamil¬ 
tonian as a matrix 

H(t) = Y^B^[K - Y(^(f)P)Y^]BY + U(t). (75) 

^Here and below all equations are expressed in atomic units 
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Here the kinetic energy operator matrix K, the potential energy operator 
matrix U and the momentum operator matrix P are diagonal, and their 
elements are written as: 


^nlmn'Um' 


(76) 

Uijki^jfk' (^) 


(77) 

Pnjkn'j'k' 

kn'^jk^nn^^jj'^kk'' 

(78) 


Hence the operations number needed for the multiplication by matrix fl75l) 
scales as the one for the transform Eq. (l6^ . In the absence of the vector 
potential the multiplication by the Hamiltonian requires only two transforms, 
namely the direct and inverse ones. When the vector potential is non-zero, 
one employs the additional couple of angular transforms Y. In the present 
case they perform the transition from the expansion in spherical waves to the 
one in plane waves (that is, to the DVR in the momentum representation) 
and vice versa. Since all the transforms are orthogonal, the Hamiltonian 
matrix preserves the original Hamiltonian hermiticity. 

Next, we take the opportunity to introduce in the kinetic energy matrix 
Eq. fl76l) the angular momentum-dependent momentum kni. One can dehne 
it either just as kni = kn, or in a more advanced fashion, namely via Eq. flbOj) . 
We will compare these two ways below. 

The TDSE has the form 


^ dt 




(79) 


In the matrix form it is written as 


dt 


= m)m- 


(80) 
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Since one may perform the fast fast multiplication by matrix of the form 
of fl75|l . the time dependent equation can be solved by means of different 
approaches, e.g. the leap-frog method or short iterative Lanczos propagator 


method [17[. However here we shall use the split-operator method 
might be represented in the form of equations in the following order: 


18|. It 



= exp[-zU(t-f r/2)r/2]'0(t); 

Clit) 

= BY7/)i(f); 

C2(t) 

= exp[—zKr/2]ci(f); 

^2{t) 

= Ytc2(t); 


= exp[i(^(t-h r/ 2 )P)r]c^ 2 (^); 

o 

CO 

= Y^3(t); 

Ci{t) 

= exp[-iKr/2]c3(t); 


= YtBtc4(t); 

+ t) 

= exp[—zU(f-f r/2)r/2]'04(f) 


This sequence of steps is equivalent to '4>{t + t) = exp[— iH(f -|- r/2)r]'j/j(t) 
with the accuracy 0{t^), so the method has a global error of O(r^). As 
the matrices K, U and P are diagonal, the exponential functions of them 
reduce to the exponential functions of complex numbers. Therefore each step 
of the method performing requires a number of operations 0{NrNgN'^) + 
0(7V,7VX) + 0{NeN^Nr log^ iV,). 

Upon the employing of the approach that is being presented, the evo¬ 
lution of the phases of free spherical waves is evaluated more precisely, the 
greater the evaluation region. It emerges to be an important advantage in 
comparison to another space approximation techniques that are commonly 
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used today (finite-difference method, finite-element method and so on). This 
is extremely helpful for the problems that require the consideration of the 
long-duration wavefunction evolution in large space regions in variable ex¬ 
ternal fields. 

It is worth mentioning that, although we are considering the TDSE solv¬ 
ing only, the reduction of a problem to the multiplication by the matrix of 
the form of d75|) might be also used in iteration methods for the stationary 
elliptic equations solving as well. 


Test on 3D time-dependent harmonic oscillator 


( 0=1 


( 0=2 




Figure 6: Convergence of the solution to the exact one for the oscillator in external time- 
dependent field: first row — the coordinate gauge; second row — the velocity gauge; left 
column — external field frequency a; = 1; right column — uj = 2. 


As a first benchmark application let us consider the problem of 3D har- 
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monic oscillator in external field, that has the analytical solution. As this 
problem possesses features somewhat opposite to those which are optimal for 
the employing of DSBT-DVR (that is, it needs only the small spatial region 
size), it proves to be the most stringent test for the method. 

The spherically-symmetric three-dimensional harmonic oscillator poten¬ 
tial is known to be 


U,{r) = 


r 


2 


2 ■ 


The time-dependent external field can be presented by means of various ways 
which are equivalent in terms of theory, but different in terms of their imple¬ 
mentation by a numerical scheme. We have accomplished the calculations 
for the external held representation both in the coordinate gauge 


V{f,t) = -q£{t)z- 
A{t) = 0. 


and the velocity gauge 


V(r, t) = 0; 
A{t) = A{t)ez. 


The pulse form was supposed to be 

A{t) = —^0 sin cat. 


and 


q£{t) 


dA 


ujAq cos ut. 
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We took the external field amplitude to be = 0.25. Computations were 
carried out for the two external field frequencies, namely a; = 1 correspond¬ 
ing to the oscillator resonant frequency, and the non-resonant u = 2. At the 
resonant frequency the amplitude of the wavepacket center position oscilla¬ 
tions (against center of coordinate) grows linearly with time, that is, higher 
and higher spherical harmonics are excited, whereas in the non-resonant case 
only small ^ harmonics are excited. Initial state function was set equal to 
the ground state function. 

We used the parameter 5{t) = |1 — ('0osc(f^, ^))| to estimate the ap¬ 

proximation error. Here wave function ‘iposc{^,t) is the analytical solution for 
the three-dimensional harmonic oscillator in a time-dependent external field, 
and is the numerical solution. We set the angular basis parameters 

No = 16 and = 1. In order to diminish the error of the split-operator 
method (which is of no current interest), the time step has been set very 
small. 

The Figure [6] shows the error 6{t) of the obtained numerical solution as a 
function of time at kni = kn for the three grids having the same step Ar = 0.2, 
but different = 12.8, 25.6, and 51.2. It is apparent that the solution 
error falls down as 0(r“^), as one should expect basing on the fact that the 
approximate transform error makes the most significant contribution to the 
overall scheme error at such scheme parameters. 

The Figure [7] presents the same as the Figure El except that we have 
employed the corrected kni from Eq. fl60|) . One can see that the rate of the 
solution convergence to the exact one depending on r^ax is the same as in 
the case kni = kn, but the error absolute value is 4 to 5 times less. 
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Figure 7: Convergence of the solution to the exact one for the oscillator in external time- 
dependent field for the corrected kni- first row — the coordinate gauge; second row — the 
velocity gauge; left column — external field frequency w = 1; right column — w = 2. 


4-3. Molecule in the pump-probe field 

Now let us turn to a problem of more physical use. As a benchmark 
example we shall consider the molecule in the held of complex-shaped 
laser pulse consisting of the short ultraviolet (XUV) pulse combined with 


the long infrared radiation (IR) pu^ 
developing pump-probe techniques 


se. This emerges as a model of rapidly 


2n[ |. The electron is emitted after being 


subjected to the XUV pumping pulse and then moves under joint action of the 
long-range Coulomb held and slowly changing IR probe pulse held. Modeling 
of this process requires computations for a long atomic time period as well 
as for the large simulation region size r^ax (in order for the electron not to 
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escape outside its boundaries). Since a large r^ax implies a small momentum 
step, an DSBT based approach emerges to be perfectly appropriate for this 
problem solving. 

First we need to estimate the accuracy that our scheme provides for 
the singular potential problems which are frequently encountered in atomic 
physics. To this end, we have computed the eigenenergies of the approximate 
Hamiltonian fl75|) for the different singular potentials. The ground state en¬ 
ergy and wavefunction have been evaluated by means of the imaginary time 
evolution method (that is to substitute t —it in Eq. fl80|) h The excited 
states have been evaluated via the imaginary time evolution method with 
the ortogonalization of the wavefunction to the lower states functions on 
each time step. 

We shall begin with the considering of the Hydrogen atom whose nucleus 
potential is known to be 



The table [T] demonstrates the convergence of the calculated energy with the 
grid step Ar decreasing at the hxed r^ax = 102.4. It is seen that calculated 
energies oi i = 0-states converge to the exact ones quadratically. Meanwhile, 
i > 0-states energies hardly depend on Ar and possess much less errors. The 
latter are caused mainly by the very SET error and decrease quadratically 
with Tjnax increasing. For £ = 0, the large error value and rather slow Ar- 
convergence result from the fact that the Coulomb wavefunctions with £ = 0 
have the hrst derivative discontinuity at r = 0 and are poorly approximated 
by the sin Fourier expansion. 

In order to enhance the convergence rate for £ = 0, one can replace the 
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Table 1: Bound states energies for H. 


n 1 

Ar 


Exact 

0.2 0.1 

0.05 

1 0 

-0.505927 -0.501575 

-0.500405 

-0.5 

2 0 

-0.125738 -0.125197 

-0.125051 

-0.125 

3 0 

-0.055774 -0.055614 

-0.055571 

-0.055555(5) 

2 1 

-0.125017 -0.125017 

-0.125017 

-0.125 

3 1 

-0.055572 -0.055572 

-0.055572 

-0.125 

3 2 

-0.055606 -0.055606 

-0.055606 

-0.055555(5) 


exact Coulomb potential with an effective potential constructed in such a 
manner that, at a given approximate kinetic energy operator Y^BlKBY, the 
approximated Hamiltonian ground state function and energy would coincide 
with the exact ones for the ground state of a hydrogen-like ion with a nucleus 
charge Z, that is, correspondingly, 

^ 3/2 ^2 

^zwo{^ = ^^exp(-Zr); Ezio = —(82) 

V TT 2 , 


For a nucleus residing in a point with the coordinates r^, such a potential is 
expressed as 


ijki ^ a) 


[YtBt[K - Eziol\BYipz{ra)]ijk 

[¥’z{ra)]ijk 


(83) 


where 


[^z{ra)]ijk = ^zioo{rijk - ra)ri\/^rArjjAcj). (84) 

However, since (pziooi^^ tends to zero exponentially at large r’s, the expres¬ 
sion fl5^ would yield the result going to infinity at large |r — ra\ due to 
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numerical errors. To avoid this, we have chosen to use the following poten¬ 
tial 


uz{rijk,ra) = f{\rijk - ra\)uz{rijk,ra) - ^ (85) 

Here f{r) is the mask function possessing the properties /(O) = 1, f{r —)■ 
cxd) = 0. In our calculations the mask function of the form 

/(r) = exp(—Zr). (86) 

was employed. In such a way, the potential (|85|1 coincides with the potential 
(|83|l when \ f—fa\ is small and with the usual Coulomb potential when it is 
large. As all the Coulomb functions at r ^ 0 have an asymptotic behavior 
~ 1 —Zr-|-(9(r^), the increasing of the accuracy of near-r = 0 approximation 
of the ground state function (pzioo(^ should lead to the increasing of the 
accuracy of the approximation of the other Coulomb functions. 

Table 2: Bound state energies for H with the effective potential in use. 




Ar 


n i 

0.2 

0.1 

0.05 

1 0 

-0.500967 

-0.500133 

-0.500017 

2 0 

-0.125125 

-0.125017 

-0.125002 

3 0 

-0.055593 

-0.055561 

-0.055556 

2 1 

-0.125016 

-0.125017 

-0.125017 

3 1 

-0.055572 

-0.055572 

-0.055572 

3 2 

-0.055606 

-0.055606 

-0.055606 


The table |2] exhibits the same as the Table [T] does, except that the po- 
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tential 


Uoi^r) = uz{f,Q) 

has been used instead of the Coulomb one. One can easily observe the 
decreasing of the differences between the calculated energies of £ = 0-states 
and the exact energies of H atom stationary states, whereas for the £ > 0- 
states these differences apparently do not increase. 

Now turn to the molecular Hydrogen ion H^. When dealing with multi- 
nuclear systems, one has to employ a potential with singularities that do not 
coincide with the coordinate origin. We shall write the approximate nuclear 
potential in the hydrogen molecule in the following way: 

Uo{r) = Mi(r, R/2) + Mi(f, -R/2). 

where R is the internuclear vector with the length R = 2 which corresponds 
to the equilibrium internuclear distance for the ground state. We have 
chosen the internuclear direction along the Oz axis orientation, R\\ez- The 



Table 3: 

Bound state 

energies for HJ. 



Ng 


State 

4 

8 

16 Exact 

lag 

-1.066449 

-1.094991 

-1.101242 -1.102634 

2,(Ju 

-0.618383 

-0.659020 

-0.666117 -0.667534 


table [T] manifests the calculated energies for ground and hrst excited 
states converge with the angular basis parameter Ng increasing at the hxed 
Ar = 0.2 and r^ax = 102.4. The “exact” energies given here were obtained 
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through the calculation via the method |l9j based upon the spheroidal co¬ 
ordinates utilizing. The error arising from the grid step is negligible in this 
case, therefore the table of convergence over the grid step is not presented 
here. 

Next let us consider the evolution of the molecular ion in the held of 
two overlapping linearly polarized laser pulses 

A{t) = Ajjv{t)njjv + AiR{t)niR. 

Here the “XUV” pulse was supposed to have the Gaussian envelope 


•^uv(t) = —v4.(7v exp ( —2 In 2 


w 


uv 


cos cot 


where wuv is the full width at half maximum . Next, the “IR” pulse was 
chosen to have a compact support and the cos^-envelope, as follows 

Ainit) = —AiRCos^[7i{t — tiR)/TiR] cosuiR{t — tiR), \t — tjR\ < tir/2, 

where tjr is the overall pulse duration, tjR is the shift of the arrival time of the 
IR-pulse center relative to that for the XUV pulse. The external held of this 


form is employed in the attosecond streaking method 


20j. The XUV-pulse 


triggers the ionization, then the detected electrons spectrum dependence on 
the time shift tjR enables to determine the IR pulse genuine form, or, in 
the case of this form being known, to obtain the time delay of the electron 
emission during the ionization process. 

The probe pulse parameters was taken to be ujr = 0.062832, Air = 0.05, 
and tir = 2Tir = 200 (which are common values in modern attosecond 
streaking experiments), and the pump pulse parameters, correspondingly, 
were urv = l-^'ol + 0.5 {Eq standing for the molecule ground state energy. 
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evaluated by means of the imaginary time evolution method) Auv = 0.25, 
and wuv = 10. Both pulses polarization were chosen to be co-directed with 
the molecular axis, fiuv = 'h/J? = ez- In all the examples referred to below 
we used the numerical scheme parameters as follows: r^ax = 409.6, time 
step r = Ar^/4, evolution beginning time to = + tin, and evolution 

termination time tfin = tir/2 = 100. 



Figure 8: (left) The probability density P = versus r at 0 = 0 for the different 

tin’s; (right) gray scale map of log^g P versus Cartesian coordinates x and z at ?/ = 0 for 
tin = 0. 


The Fig. [H] shows the probability density P{f,tfin) = |4'(r, that 

the electron is at r. The calculations were performed for the scheme param¬ 
eters Ng = 16 and Ar = 0.2. Due to the stationary phase approximation, 
for the time tfin S> wuv and for r S> 1, the relation P(f,tfin) ~ (rif/tfin) 
holds, where <j{k) is the differential cross section of electrons emission de¬ 
pending on momentum. On the left panel of the Fig. [8], the peak near 
r = 0 corresponds to the wavefunction of the ground and other stationary 
states, whereas the peak centered in the vicinity of r = 100 emerges due to 
the one-photon ionization, and and the rest large r peaks are caused by the 
multiphoton processes. This is apparently conhrmed by the right panel of 
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the Fig. El where the r = 100 enhanced probability ring has one node de¬ 
pending on 6, the circle of larger radius has two nodes corresponding to the 
dipole and quadrupole distributions arising from the absorption of one or two 
photons correspondingly. The left panel of the Figure [8] also demonstrates 
the probability density dependence on the IR pulse phase at the moment of 
the XUV pulse arrival. The theory predicts the probe pulse action causing 
the electron momentum shift equal to the magnitude of the IR pulse at the 
moment of the electron emission from the molecule (which roughly coincides 
with the moment of the XUV pulse arrival). This is exactly what is observed 
on the right panel of the Figure El 
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Figure 9: (left) The pairwise differences A = \P 2 Ar,Ng — PAr,Ne\i"^ at the fixed Ng = 8; 
(right) The pairwise differences A = \PAr,Ne /2 — PAr,Ne\r'^ at the fixed Ar = 0.2. 


Besides, we have examined the probability density P convergence rate 
depending on the step Ar and on the angular basis size, Ng. The left 
panel of the |9] presents the pairwise differences A(r) = \P 2 Ar,Ng{,rez,tfin) — 
PAr,Ng{rez,tfin)\r^ for the probability densities PAr,Ng{,r,tfin) evaluated on 
the three grids having the steps Ar = 0.4, 0.2 and 0.1 at the hxed Ng = 8. 
For the sake of comparison PAr,Ng{rez,tfin)r‘^ for Ar = 0.1 and Ng = 8 is 
plotted on the same hgure. It is apparent that even on the coarsest grid 
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with Ar = 0.4 the error is of the order of 1%; this value is quite small in 
terms of experimental accuracy which is common in the held in question. 

Upon halving the step size, the error drops down by 1-2 orders of magnitude. 

However, the error in a particular point decreases non-uniformly, actually as 
expected due to the global basis functions using. 

The right panel of the[9]displays the pairwise differences A(r) = \PAr,Ne/ 2 {'rez,tfin) — 
PAr,Ng{'rez,tfin)\r‘^ for the three different angular bases with Nq = 4, 8, and 
16 at the hxed step Ar = 0.2, as well as PAr,Ng{'i"^Zi'tfin)'i"‘^ for Ar = 0.2 and 
Ng = 16. One can see that the error due to the angular basis small size is 
much larger than that due to the radius step. This is related to the molecular 
potential non-centrality. For Ng = 4 the error has magnitude about 25% (in 
the vicinity of maxima), whereas upon the basis size increasing up to = 8 
the error drops down to 6%. Therefore Ng = 16 has been chosen for the main 
part of our calculations. 

5. Conclusion 

We have developed the algorithm for the DSBT that possesses the advan¬ 
tages of orthogonality, performing fastness and uniform grid. Our approach is 
based upon the SET factorization into the two subsequent orthogonal trans¬ 
forms, namely the fast Fourier transform (requiring the operations number 
0(Alog2 A)) and the orthogonal transform founded on the discrete orthog¬ 
onal Legendre polynomials (requiring the operations number 0{iN)). Our 
discrete transform converges to the exact SET as the square of the momen¬ 
tum grid step. 

Besides, basing on DSBT and DVR, we have also elaborated the 3D 
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TDSE solving method (DSBT-DVR). The examination of the DSBT-DVR 
algorithm has demonstrated its efficiency for the purposes of solving of time- 
dependent problems in atomic and molecular physics. An DSBT based ap¬ 
proach allows to evaluate the free spherical wave functions evolution the more 
accurately, the more is the spatial region size. It appears to be an advantage 
in comparison to another methods applied in this held. This is especially 
helpful for problems like the modelling of the attosecond streaking approach 
and other pump-probe techniques, since they require the computation of the 
wavefunction evolution under the joint action of long-lasting pulses and the 
weak Coulomb held on large spatial regions. Another important preference 
of the method proposed is the fast convergence over grid step when applied 
to the problems with smooth (or artihcially smoothed) potentials. 

It should be noted that the current DSBT-DVR version does not make 
any use of another helpful DBBT feature, namely the DSBT capability to 
be employed for the aim of the evaluation of multi-center integrals 6j. The 
leveraging of this capability for the solving of both SSE and TDSE for the 
multielectron molecules is expected to be the matter of our future work. 
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Appendix A. Proof of transform orthogonality 


Let us begin with the demonstration of the Eq. fl32p validity. Consider the 
sum 


N 

a, = 5^V[P,](z,iV-lK (A.l) 

i=0 

It may be transformed as follows 

N N 

= ^ -hN- 1)P 

7V-1 

= ^ - 1)** + N - 1)A^' 

i=0 

N 

- Y, PS -l,N- 1)P - N - 1)0* 

i=l 

According to (El]), if s < £ the sums in the last string equal zero, hence 

a(s) = Pii-1, A - 1) [(-l)'A* - 0*] . (A.2) 

Here we also used the DLOP parity property 

Pi{t, N) = {-ifPeiN -t,N). (A.3) 

Now let us take an arbitrary discrete polynomial 

pit) = YCsP-^ (A.4) 

s=0 

and consider the sum 

N M 

a = V[Pi]{t, N - l)p{t) = Y Csas (A.5) 

i=0 s=0 
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By using Eq. flA.2p we obtain 


<7 = n(-i,]v-i)[(-i)'p(iv)-p(o)] 


(A.6) 


Next, one can construct the weighted sum 


N 


J2nPi]ihN-l)p{z)wi{N) = 


a - ^ V[P,](0, N - l)p(O) - ^V[P,](iV, N - l)p(iV) 


Making use of the Eqs. flA.bi [35l IA.3p and the normalization condition Pi{0, N— 
1) = 1 yields 


N 




1 + Pij-^, N - 1) 


[(-IYp(N) - piO)] . 


(A.7) 


After the division of both sides of this equation by [1 + P£(—1, N — l)]/2, we 
arrive to Eq. fl3^ . 

Now let us prove Eq. fHT]) . As this equation is symmetric with respect to 
the exchange of indices n and m, for dehniteness we shall assume n > m. 
Since Lmi = 0 for / > m, one can write 



(A.8) 


For sake of the notation simplicity, from now on we designate 



(A.9) 


So, according to Eq. flA.Sp . Eq. fHTj) holds true, when 




mm 


L 


(A.IO) 
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The substitution of L elements definition from Eq. fl39p yields 

Km = ^ Piin - I, 2u)p;(m - I, 2m) (^1 “ y “ • (A-ll) 

By change of the summation index i = m — I we can rewrite flA.lip as 


^n.m. 


Pt X A A \ 

P'^{n -m + i, 2n)P'{i, 2m) M - ^ - y j . (A.12) 

i=0 ^ ' 


(A. 13) 


Due to Eq. flA.3p . DDLOP have the parity property 

The sum in Eq. (lA.12p might be split into the two sums as 

Km = ^ Pejn -m + i, 2n)P^{i, 2m) - y j 

Here we made use of the fact that P^'(m,2m) = 0 when = 1 at even 
Next we apply the parity property to both DDLOPs in the summand of the 
second sum in Eq. flA.14l) and make the summation index change i —)■ 2m — i 


ra—1 


PiK — m + i, 2n)P^{i, 2m) ( 1 — 


Si, 


iO 


i=0 

m—1 




” ^ ^ / A. \ 

Y^ P'lK + rn — i, 2n)P'^{2m — i, 2m) ( 1-^ j 

i=0 ^ 

2m / ^ 

Y + m-i, 2n)P[{i, 2m) ( 1 - 


2=m+l 


The latter sum then might be combined with the (remained unchanged) first 
sum in Eq. flA.14p to get 


1 


2m 


Km ^ 3 XI ^ 2n)P/(i, 2m)w,{2m). 


(A. 16) 


i=0 
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Ks P'(^{n — m + i, 2n) is the polynomial of the order £ — 1, we can apply Eq. fl32p 
to obtain 


2m 


^ P^{n — m + i, 2n)P^{i, 2m)wi{2m) 


i=0 


= (—+ m, 2n) — Pi{n — m, 2n). 


Finally, after using the parity property Eq. flA.13p . the result becomes 


^nm Tn^2lPj. 


(A.16) 


Since Lnm = —Pe{n — m,2n) for n > m > hqi > 0, we have thus proved 
Eq. flA.lOD and therefore Eq. fHB . 
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